/* -*- mode: C -*-  */
/* 
   IGraph library.
   Copyright (C) 2009-2012  Gabor Csardi <csardi.gabor@gmail.com>
   334 Harvard street, Cambridge, MA 02139 USA
   
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.
   
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
   
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 
   02110-1301 USA

*/

#include "config.h"

#include "cs/cs.h"

#include "igraph_sparsemat.h"
#include "igraph_error.h"
#include "igraph_interface.h"
#include "igraph_constructors.h"
#include "igraph_memory.h"
#include "igraph_vector_ptr.h"
#include "igraph_attributes.h"

#include <string.h>

/** 
 * \section about_sparsemat About sparse matrices 
 * 
 * <para>
 * The <code>igraph_sparsemat_t</code> data type stores sparse matrices, 
 * i.e. matrices in which the majority of the elements are zero.
 * </para>
 *
 * <para>The data type is essentially a wrapper to some of the
 * functions in the CXSparse library, by Tim Davis, see
 * http://www.cise.ufl.edu/research/sparse/CXSparse/
 * </para>
 * 
 * <para>
 * Matrices can be stored in two formats: triplet and
 * column-compressed. The triplet format is intended for sparse matrix
 * initialization, as it is easy to add new (non-zero) elements to
 * it. Most of the computations are done on sparse matrices in
 * column-compressed format, after the user has converted the triplet
 * matrix to column-compressed, via \ref igraph_sparsemat_compress().
 * </para>
 * 
 * <para>
 * Both formats are dynamic, in the sense that new elements can be
 * added to them, possibly resulting the allocation of more memory.
 * </para>
 * 
 * <para>
 * Row and column indices follow the C convention and are zero-based.
 * </para>
 * 
 * <para>
 * \example examples/simple/igraph_sparsemat.c
 * \example examples/simple/igraph_sparsemat2.c
 * \example examples/simple/igraph_sparsemat3.c
 * \example examples/simple/igraph_sparsemat4.c
 * \example examples/simple/igraph_sparsemat5.c
 * \example examples/simple/igraph_sparsemat6.c
 * \example examples/simple/igraph_sparsemat7.c
 * \example examples/simple/igraph_sparsemat8.c
 * </para>
 */

/**
 * \function igraph_sparsemat_init
 * Initialize a sparse matrix, in triplet format
 * 
 * This is the most common way to create a sparse matrix, together
 * with the \ref igraph_sparsemat_entry() function, which can be used to
 * add the non-zero elements one by one. Once done, the user can call
 * \ref igraph_sparsemat_compress() to convert the matrix to
 * column-compressed, to allow computations with it.
 * 
 * </para><para>The user must call \ref igraph_sparsemat_destroy() on
 * the matrix to deallocate the memory, once the matrix is no more
 * needed.
 * \param A Pointer to a not yet initialized sparse matrix. 
 * \param rows The number of rows in the matrix.
 * \param cols The number of columns.
 * \param nzmax The maximum number of non-zero elements in the
 *    matrix. It is not compulsory to get this right, but it is
 *    useful for the allocation of the proper amount of memory.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_init(igraph_sparsemat_t *A, int rows, int cols, int nzmax) {

  if (rows < 0) { 
    IGRAPH_ERROR("Negative number of rows", IGRAPH_EINVAL);
  }
  if (cols < 0) {
    IGRAPH_ERROR("Negative number of columns", IGRAPH_EINVAL);
  }
  
  A->cs=cs_spalloc( rows, cols, nzmax, /*values=*/ 1, 
		  /*triplet=*/ 1);
  if (!A->cs) {
    IGRAPH_ERROR("Cannot allocate memory for sparse matrix", IGRAPH_ENOMEM);
  }

  return 0;
}

/** 
 * \function igraph_sparsemat_copy
 * Copy a sparse matrix
 * 
 * Create a sparse matrix object, by copying another one. The source
 * matrix can be either in triplet or column-compressed format.
 * 
 * </para><para>
 * Exactly the same amount of memory will be allocated to the
 * copy matrix, as it is currently for the original one.
 * \param to Pointer to an uninitialized sparse matrix, the copy will
 *    be created here.
 * \param from The sparse matrix to copy.
 * \return Error code.
 *
 * Time complexity: O(n+nzmax), the number of columns plus the maximum
 * number of non-zero elements.
 */

int igraph_sparsemat_copy(igraph_sparsemat_t *to, 
			  const igraph_sparsemat_t *from) {

  int ne=from->cs->nz == -1 ? from->cs->n+1 : from->cs->nzmax;

  to->cs = cs_spalloc(from->cs->m, from->cs->n, from->cs->nzmax, 
		      /*values=*/ 1, 
		      /*triplet=*/ igraph_sparsemat_is_triplet(from));

  to->cs->nzmax = from->cs->nzmax;
  to->cs->m     = from->cs->m;
  to->cs->n     = from->cs->n;
  to->cs->nz    = from->cs->nz;

  memcpy(to->cs->p, from->cs->p, sizeof(int) * (size_t) ne);
  memcpy(to->cs->i, from->cs->i, sizeof(int) * (size_t) (from->cs->nzmax));
  memcpy(to->cs->x, from->cs->x, sizeof(double) * (size_t) (from->cs->nzmax));

  return 0;
}

/**
 * \function igraph_sparsemat_destroy
 * Deallocate memory used by a sparse matrix
 * 
 * One destroyed, the sparse matrix must be initialized again, before
 * calling any other operation on it.
 * \param A The sparse matrix to destroy.
 * 
 * Time complexity: O(1).
 */

void igraph_sparsemat_destroy(igraph_sparsemat_t *A) {
  cs_spfree(A->cs);
}

/**
 * \function igraph_sparsemat_realloc
 * Allocate more (or less) memory for a sparse matrix
 * 
 * Sparse matrices automatically allocate more memory, as needed. To
 * control memory allocation, the user can call this function, to
 * allocate memory for a given number of non-zero elements.
 * \param A The sparse matrix, it can be in triplet or
 *    column-compressed format.
 * \param nzmax The new maximum number of non-zero elements.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_realloc(igraph_sparsemat_t *A, int nzmax) {
  return !cs_sprealloc(A->cs, nzmax);
}

/**
 * \function igraph_sparsemat_nrow
 * Number of rows
 * 
 * \param A The input matrix, in triplet or column-compressed format.
 * \return The number of rows in the \p A matrix.
 * 
 * Time complexity: O(1).
 */

long int igraph_sparsemat_nrow(const igraph_sparsemat_t *A) {
  return A->cs->m;
}

/**
 * \function igraph_sparsemat_ncol
 * Number of columns.
 * 
 * \param A The input matrix, in triplet or column-compressed format.
 * \return The number of columns in the \p A matrix.
 * 
 * Time complexity: O(1).
 */

long int igraph_sparsemat_ncol(const igraph_sparsemat_t *A) {
  return A->cs->n;
}

/**
 * \function igraph_sparsemat_type
 * Type of a sparse matrix (triplet or column-compressed)
 * 
 * Gives whether a sparse matrix is stored in the triplet format or in
 * column-compressed format.
 * \param A The input matrix.
 * \return Either \c IGRAPH_SPARSEMAT_CC or \c
 * IGRAPH_SPARSEMAT_TRIPLET.
 * 
 * Time complexity: O(1).
 */

igraph_sparsemat_type_t igraph_sparsemat_type(const igraph_sparsemat_t *A) {
  return A->cs->nz < 0 ? IGRAPH_SPARSEMAT_CC : IGRAPH_SPARSEMAT_TRIPLET;
}

/**
 * \function igraph_sparsemat_is_triplet
 * Is this sparse matrix in triplet format?
 * 
 * Decides whether a sparse matrix is in triplet format. 
 * \param A The input matrix.
 * \return One if the input matrix is in triplet format, zero
 * otherwise.
 * 
 * Time complexity: O(1).
 */

igraph_bool_t igraph_sparsemat_is_triplet(const igraph_sparsemat_t *A) {
  return A->cs->nz >= 0;
}

/**
 * \function igraph_sparsemat_is_cc
 * Is this sparse matrix in column-compressed format?
 * 
 * Decides whether a sparse matrix is in column-compressed format. 
 * \param A The input matrix.
 * \return One if the input matrix is in column-compressed format, zero
 * otherwise.
 * 
 * Time complexity: O(1).
 */

igraph_bool_t igraph_sparsemat_is_cc(const igraph_sparsemat_t *A) {
  return A->cs->nz < 0;
}

/**
 * \function igraph_sparsemat_permute
 * Permute the rows and columns of a sparse matrix
 * 
 * \param A The input matrix, it must be in column-compressed format.
 * \param p Integer vector, giving the permutation of the rows.
 * \param q Integer vector, the permutation of the columns.
 * \param res Pointer to an uninitialized sparse matrix, the result is
 *   stored here.
 * \return Error code.
 * 
 * Time complexity: O(m+n+nz), the number of rows plus the number of
 * columns plus the number of non-zero elements in the matrix.
 */

int igraph_sparsemat_permute(const igraph_sparsemat_t *A,
			     const igraph_vector_int_t *p, 
			     const igraph_vector_int_t *q,
			     igraph_sparsemat_t *res) {

  long int nrow=A->cs->m, ncol=A->cs->n;
  igraph_vector_int_t pinv;
  long int i;

  if (nrow != igraph_vector_int_size(p)) {
    IGRAPH_ERROR("Invalid row permutation length", IGRAPH_FAILURE);
  }
  if (ncol != igraph_vector_int_size(q)) {
    IGRAPH_ERROR("Invalid column permutation length", IGRAPH_FAILURE);
  }

  /* We invert the permutation by hand */
  IGRAPH_CHECK(igraph_vector_int_init(&pinv, nrow));
  IGRAPH_FINALLY(igraph_vector_int_destroy, &pinv);
  for (i=0; i<nrow; i++) {
    VECTOR(pinv)[ VECTOR(*p)[i] ] = (int) i;
  }
  
  /* And call the permutation routine */
  if (! (res->cs = cs_permute(A->cs, VECTOR(pinv), VECTOR(*q), /*values=*/ 1))) {
    IGRAPH_ERROR("Cannot index sparse matrix", IGRAPH_FAILURE);
  }
  
  igraph_vector_int_destroy(&pinv);
  IGRAPH_FINALLY_CLEAN(1);
  
  return 0;
}

int igraph_i_sparsemat_index_rows(const igraph_sparsemat_t *A,
				  const igraph_vector_int_t *p,
				  igraph_sparsemat_t *res,
				  igraph_real_t *constres) {

  igraph_sparsemat_t II, II2;
  long int nrow=A->cs->m;
  long int idx_rows=igraph_vector_int_size(p);
  long int k;

  /* Create index matrix */
  IGRAPH_CHECK(igraph_sparsemat_init(&II2, (int) idx_rows, (int) nrow, 
				     (int) idx_rows));
  IGRAPH_FINALLY(igraph_sparsemat_destroy, &II2);
  for (k=0; k<idx_rows; k++) {
    igraph_sparsemat_entry(&II2, (int) k, VECTOR(*p)[k], 1.0);
  }
  IGRAPH_CHECK(igraph_sparsemat_compress(&II2, &II));
  igraph_sparsemat_destroy(&II2);
  IGRAPH_FINALLY_CLEAN(1);
  IGRAPH_FINALLY(igraph_sparsemat_destroy, &II);
  
  /* Multiply */
  IGRAPH_CHECK(igraph_sparsemat_multiply(&II, A, res));
  igraph_sparsemat_destroy(&II);
  IGRAPH_FINALLY_CLEAN(1);
  
  if (constres) {
    if (res->cs->p[1] != 0) {
      *constres = res->cs->x[0];
    } else {
      *constres = 0.0;
    }
  }

  return 0;
}

int igraph_i_sparsemat_index_cols(const igraph_sparsemat_t *A,
				  const igraph_vector_int_t *q,
				  igraph_sparsemat_t *res,
				  igraph_real_t *constres) {

  igraph_sparsemat_t JJ, JJ2;
  long int ncol=A->cs->n;
  long int idx_cols=igraph_vector_int_size(q);
  long int k;
  
  /* Create index matrix */
  IGRAPH_CHECK(igraph_sparsemat_init(&JJ2, (int) ncol, (int) idx_cols, 
				     (int) idx_cols));
  IGRAPH_FINALLY(igraph_sparsemat_destroy, &JJ2);
  for (k=0; k<idx_cols; k++) {
    igraph_sparsemat_entry(&JJ2, VECTOR(*q)[k], (int) k, 1.0);
  }
  IGRAPH_CHECK(igraph_sparsemat_compress(&JJ2, &JJ));
  igraph_sparsemat_destroy(&JJ2);
  IGRAPH_FINALLY_CLEAN(1);
  IGRAPH_FINALLY(igraph_sparsemat_destroy, &JJ);

  /* Multiply */
  IGRAPH_CHECK(igraph_sparsemat_multiply(A, &JJ, res));
  igraph_sparsemat_destroy(&JJ);
  IGRAPH_FINALLY_CLEAN(1);

  if (constres) {
    if (res->cs->p [1] != 0) {
      *constres = res->cs->x [0];
    } else {
      *constres = 0.0;
    }
  }

  return 0;
}

/**
 * \function igraph_sparsemat_index
 * Index a sparse matrix, extract a submatrix, or a single element
 * 
 * This function serves two purposes. First, it can extract
 * submatrices from a sparse matrix. Second, as a special case, it can
 * extract a single element from a sparse matrix.
 * \param A The input matrix, it must be in column-compressed format.
 * \param p An integer vector, or a null pointer. The selected row
 *    index or indices. A null pointer selects all rows.
 * \param q An integer vector, or a null pointer. The selected column
 *    index or indices. A null pointer selects all columns.
 * \param res Pointer to an uninitialized sparse matrix, or a null
 *    pointer. If not a null pointer, then the selected submatrix is
 *    stored here.
 * \param constres Pointer to a real variable or a null pointer. If
 *    not a null pointer, then the first non-zero element in the
 *    selected submatrix is stored here, if there is one. Otherwise
 *    zero is stored here. This behavior is handy if one
 *    wants to select a single entry from the matrix.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_index(const igraph_sparsemat_t *A,
			   const igraph_vector_int_t *p,
			   const igraph_vector_int_t *q,
			   igraph_sparsemat_t *res,
			   igraph_real_t *constres) {

  igraph_sparsemat_t II, JJ, II2, JJ2, tmp;
  long int nrow=A->cs->m;
  long int ncol=A->cs->n;
  long int idx_rows= p ? igraph_vector_int_size(p) : -1;
  long int idx_cols= q ? igraph_vector_int_size(q) : -1;
  long int k;

  igraph_sparsemat_t *myres=res, mres;

  if (!p && !q) {
    IGRAPH_ERROR("No index vectors", IGRAPH_EINVAL);
  }

  if (!res && (idx_rows != 1 || idx_cols != 1)) {
    IGRAPH_ERROR("Sparse matrix indexing: must give `res' if not a "
		 "single element is selected", IGRAPH_EINVAL);
  }

  if (!q) {
    return igraph_i_sparsemat_index_rows(A, p, res, constres);
  }
  if (!p) {
    return igraph_i_sparsemat_index_cols(A, q, res, constres);
  }

  if (!res) {
    myres=&mres;
  }

  /* Create first index matrix */
  IGRAPH_CHECK(igraph_sparsemat_init(&II2, (int) idx_rows, (int) nrow, 
				     (int) idx_rows));
  IGRAPH_FINALLY(igraph_sparsemat_destroy, &II2);
  for (k=0; k<idx_rows; k++) {
    igraph_sparsemat_entry(&II2, (int) k, VECTOR(*p)[k], 1.0);
  }
  IGRAPH_CHECK(igraph_sparsemat_compress(&II2, &II));
  igraph_sparsemat_destroy(&II2);
  IGRAPH_FINALLY_CLEAN(1);
  IGRAPH_FINALLY(igraph_sparsemat_destroy, &II);

  /* Create second index matrix */
  IGRAPH_CHECK(igraph_sparsemat_init(&JJ2, (int) ncol, (int) idx_cols, 
				     (int) idx_cols));
  IGRAPH_FINALLY(igraph_sparsemat_destroy, &JJ2);
  for (k=0; k<idx_cols; k++) {
    igraph_sparsemat_entry(&JJ2, VECTOR(*q)[k], (int) k, 1.0);
  }
  IGRAPH_CHECK(igraph_sparsemat_compress(&JJ2, &JJ));
  igraph_sparsemat_destroy(&JJ2);
  IGRAPH_FINALLY_CLEAN(1);
  IGRAPH_FINALLY(igraph_sparsemat_destroy, &JJ);
  
  /* Multiply */
  IGRAPH_CHECK(igraph_sparsemat_multiply(&II, A, &tmp));
  igraph_sparsemat_destroy(&II);
  IGRAPH_FINALLY_CLEAN(1);
  IGRAPH_FINALLY(igraph_sparsemat_destroy, &tmp);
  IGRAPH_CHECK(igraph_sparsemat_multiply(&tmp, &JJ, myres));
  igraph_sparsemat_destroy(&tmp);
  igraph_sparsemat_destroy(&JJ);
  IGRAPH_FINALLY_CLEAN(2);

  if (constres) {
    if (myres->cs->p [1] != 0) {
      *constres = myres->cs->x [0];
    } else {
      *constres = 0.0;
    }
  }

  if (!res) {
    igraph_sparsemat_destroy(myres);
  }

  return 0;
}

/**
 * \function igraph_sparsemat_entry
 * Add an element to a sparse matrix
 * 
 * This function can be used to add the entries to a sparse matrix,
 * after initializing it with \ref igraph_sparsemat_init().
 * \param A The input matrix, it must be in triplet format.
 * \param row The row index of the entry to add.
 * \param col The column index of the entry to add.
 * \param elem The value of the entry.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_entry(igraph_sparsemat_t *A, int row, int col, 
			   igraph_real_t elem) {
  
  if (!cs_entry(A->cs, row, col, elem)) {
    IGRAPH_ERROR("Cannot add entry to sparse matrix", 
		 IGRAPH_FAILURE);
  }

  return 0;
}

/**
 * \function igraph_sparsemat_compress
 * Compress a sparse matrix, i.e. convert it to column-compress format
 * 
 * Almost all sparse matrix operations require that the matrix is in
 * column-compressed format.
 * \param A The input matrix, it must be in triplet format.
 * \param res Pointer to an uninitialized sparse matrix object, the
 *    compressed version of \p A is stored here.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_compress(const igraph_sparsemat_t *A, 
			      igraph_sparsemat_t *res) {

  if (! (res->cs=cs_compress(A->cs)) ) {
    IGRAPH_ERROR("Cannot compress sparse matrix", IGRAPH_FAILURE);
  }

  return 0;
}

/**
 * \function igraph_sparsemat_transpose
 * Transpose a sparse matrix
 * 
 * \param A The input matrix, column-compressed or triple format.
 * \param res Pointer to an uninitialized sparse matrix, the result is
 *    stored here.
 * \param values If this is non-zero, the matrix transpose is
 *    calculated the normal way. If it is zero, then only the pattern
 *    of the input matrix is stored in the result, the values are not.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_transpose(const igraph_sparsemat_t *A, 
			       igraph_sparsemat_t *res, 
			       int values) {

  if (A->cs->nz < 0) {
    /* column-compressed */
    if (! (res->cs=cs_transpose(A->cs, values)) ) {
      IGRAPH_ERROR("Cannot transpose sparse matrix", IGRAPH_FAILURE); 
    }
  } else {
    /* triplets */
    int *tmp;
    IGRAPH_CHECK(igraph_sparsemat_copy(res, A));
    tmp = res->cs->p;
    res->cs->p = res->cs->i;
    res->cs->i = tmp;
  }
  return 0;
}

igraph_bool_t 
igraph_i_sparsemat_is_symmetric_cc(const igraph_sparsemat_t *A) {
  igraph_sparsemat_t t, tt;
  igraph_bool_t res;
  int nz;

  IGRAPH_CHECK(igraph_sparsemat_transpose(A, &t, /*values=*/ 1));
  IGRAPH_FINALLY(igraph_sparsemat_destroy, &t);
  IGRAPH_CHECK(igraph_sparsemat_dupl(&t));
  IGRAPH_CHECK(igraph_sparsemat_transpose(&t, &tt, /*values=*/ 1));
  igraph_sparsemat_destroy(&t);
  IGRAPH_FINALLY_CLEAN(1);
  IGRAPH_FINALLY(igraph_sparsemat_destroy, &tt);
  IGRAPH_CHECK(igraph_sparsemat_transpose(&tt, &t, /*values=*/ 1));
  IGRAPH_FINALLY(igraph_sparsemat_destroy, &t);

  nz=t.cs->p[t.cs->n];
  res = memcmp(t.cs->i, tt.cs->i, sizeof(int) * (size_t) nz) == 0;
  res = res && memcmp(t.cs->p, tt.cs->p, sizeof(int) *
		      (size_t)(t.cs->n+1)) == 0;
  res = res && memcmp(t.cs->x, tt.cs->x, sizeof(igraph_real_t) * (size_t)nz)==0;

  igraph_sparsemat_destroy(&t);
  igraph_sparsemat_destroy(&tt);
  IGRAPH_FINALLY_CLEAN(2);

  return res;
}

igraph_bool_t 
igraph_i_sparsemat_is_symmetric_triplet(const igraph_sparsemat_t *A) {
  igraph_sparsemat_t tmp;
  igraph_bool_t res;
  IGRAPH_CHECK(igraph_sparsemat_compress(A, &tmp));
  IGRAPH_FINALLY(igraph_sparsemat_destroy, &tmp);
  res=igraph_i_sparsemat_is_symmetric_cc(&tmp);
  igraph_sparsemat_destroy(&tmp);
  IGRAPH_FINALLY_CLEAN(1);
  return res;
}

igraph_bool_t igraph_sparsemat_is_symmetric(const igraph_sparsemat_t *A) {

  if (A->cs->m != A->cs->n) { return 0; }

  if (A->cs->nz < 0) {
    return igraph_i_sparsemat_is_symmetric_cc(A);
  } else {
    return igraph_i_sparsemat_is_symmetric_triplet(A);
  }
}

/**
 * \function igraph_sparsemat_dupl
 * Remove duplicate elements from a sparse matrix
 * 
 * It is possible that a column-compressed sparse matrix stores a
 * single matrix entry in multiple pieces. The entry is then the sum
 * of all its pieces. (Some functions create matrices like this.) This
 * function eliminates the multiple pieces.
 * \param A The input matrix, in column-compressed format.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_dupl(igraph_sparsemat_t *A) {

  if (!cs_dupl(A->cs)) {
    IGRAPH_ERROR("Cannot remove duplicates from sparse matrix", 
		 IGRAPH_FAILURE);
  }
  
  return 0;
}

/**
 * \function igraph_sparsemat_fkeep
 * Filter the elements of a sparse matrix
 * 
 * This function can be used to filter the (non-zero) elements of a
 * sparse matrix. For all entries, it calls the supplied function and
 * depending on the return values either keeps, or deleted the element
 * from the matrix.
 * \param A The input matrix, in column-compressed format.
 * \param fkeep The filter function. It must take four arguments: the
 *    first is an \c int, the row index of the entry, the second is 
 *    another \c int, the column index. The third is \c igraph_real_t,
 *    the value of the entry. The fourth element is a \c void pointer,
 *    the \p other argument is passed here. The function must return
 *    an \c int. If this is zero, then the entry is deleted, otherwise
 *    it is kept.
 * \param other A \c void pointer that is passed to the filtering
 * function.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_fkeep(igraph_sparsemat_t *A, 
			   int (*fkeep)(int, int, igraph_real_t, void*),
			   void *other) {
  
  if (!cs_fkeep(A->cs, fkeep, other)) {
    IGRAPH_ERROR("Cannot filter sparse matrix", IGRAPH_FAILURE);
  }

  return 0;
}

/**
 * \function igraph_sparsemat_dropzeros
 * Drop the zero elements from a sparse matrix
 * 
 * As a result of matrix operations, some of the entries in a sparse
 * matrix might be zero. This function removes these entries.
 * \param A The input matrix, it must be in column-compressed format.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_dropzeros(igraph_sparsemat_t *A) {

  if (!cs_dropzeros(A->cs)) {
    IGRAPH_ERROR("Cannot drop zeros from sparse matrix", IGRAPH_FAILURE);
  }

  return 0;
}

/**
 * \function igraph_sparsemat_droptol
 * Drop the almost zero elements of a sparse matrix
 * 
 * This function is similar to \ref igraph_sparsemat_dropzeros(), but it
 * also drops entries that are closer to zero than the given tolerance
 * threshold.
 * \param A The input matrix, it must be in column-compressed format.
 * \param tol Real number, giving the tolerance threshold.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_droptol(igraph_sparsemat_t *A, igraph_real_t tol) {
  
  if (!cs_droptol(A->cs, tol)) {
    IGRAPH_ERROR("Cannot drop (almost) zeros from sparse matrix", 
		 IGRAPH_FAILURE);
  }
  
  return 0;
}

/**
 * \function igraph_sparsemat_multiply
 * Matrix multiplication
 * 
 * Multiplies two sparse matrices.
 * \param A The first input matrix (left hand side), in
 *   column-compressed format.
 * \param B The second input matrix (right hand side), in
 *   column-compressed format.
 * \param res Pointer to an uninitialized sparse matrix, the result is
 *   stored here.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_multiply(const igraph_sparsemat_t *A,
			      const igraph_sparsemat_t *B,
			      igraph_sparsemat_t *res) {

  if (! (res->cs=cs_multiply(A->cs, B->cs))) {
    IGRAPH_ERROR("Cannot multiply matrices", IGRAPH_FAILURE);
  }

  return 0;
}

/**
 * \function igraph_sparsemat_add
 * Sum of two sparse matrices
 * 
 * \param A The first input matrix, in column-compressed format.
 * \param B The second input matrix, in column-compressed format.
 * \param alpha Real scalar, \p A is multiplied by \p alpha before the
 *    addition. 
 * \param beta Real scalar, \p B is multiplied by \p beta before the
 *    addition. 
 * \param res Pointer to an uninitialized sparse matrix, the result
 *    is stored here.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_add(const igraph_sparsemat_t *A, 
			 const igraph_sparsemat_t *B,
			 igraph_real_t alpha,
			 igraph_real_t beta,
			 igraph_sparsemat_t *res) {

  if (! (res->cs=cs_add(A->cs, B->cs, alpha, beta))) {
    IGRAPH_ERROR("Cannot add matrices", IGRAPH_FAILURE);
  }
  
  return 0;
}

/**
 * \function igraph_sparsemat_gaxpy
 * Matrix-vector product, added to another vector.
 * 
 * \param A The input matrix, in column-compressed format.
 * \param x The input vector, its size must match the number of
 *    columns in \p A.
 * \param res This vector is added to the matrix-vector product 
 *    and it is overwritten by the result.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_gaxpy(const igraph_sparsemat_t *A,
			   const igraph_vector_t *x,
			   igraph_vector_t *res) {

  if (A->cs->n != igraph_vector_size(x) || 
      A->cs->m != igraph_vector_size(res)) {
    IGRAPH_ERROR("Invalid matrix/vector size for multiplication",
		 IGRAPH_EINVAL);
  }
  
  if (! (cs_gaxpy(A->cs, VECTOR(*x), VECTOR(*res)))) {
    IGRAPH_ERROR("Cannot perform sparse matrix vector multiplication",
		 IGRAPH_FAILURE);
  }
  
  return 0;
}

/**
 * \function igraph_sparsemat_lsolve
 * Solve a lower-triangular linear system
 * 
 * Solve the Lx=b linear equation system, where the L coefficient
 * matrix is square and lower-triangular, with a zero-free diagonal.
 * \param L The input matrix, in column-compressed format.
 * \param b The right hand side of the linear system.
 * \param res An initialized vector, the result is stored here.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_lsolve(const igraph_sparsemat_t *L,
			    const igraph_vector_t *b,
			    igraph_vector_t *res) {

  if (L->cs->m != L->cs->n) {
    IGRAPH_ERROR("Cannot perform lower triangular solve", IGRAPH_NONSQUARE);
  }

  if (res != b) {
    IGRAPH_CHECK(igraph_vector_update(res, b));
  }

  if (! cs_lsolve(L->cs, VECTOR(*res))) {
    IGRAPH_ERROR("Cannot perform lower triangular solve", IGRAPH_FAILURE);
  }
  
  return 0;
}

/**
 * \function igraph_sparsemat_ltsolve
 * Solve an upper-triangular linear system
 * 
 * Solve the L'x=b linear equation system, where the L
 * matrix is square and lower-triangular, with a zero-free diagonal.
 * \param L The input matrix, in column-compressed format.
 * \param b The right hand side of the linear system.
 * \param res An initialized vector, the result is stored here.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_ltsolve(const igraph_sparsemat_t *L,
			     const igraph_vector_t *b,
			     igraph_vector_t *res) {
  
  if (L->cs->m != L->cs->n) {
    IGRAPH_ERROR("Cannot perform transposed lower triangular solve",
		 IGRAPH_NONSQUARE);
  }
  
  if (res != b) {
    IGRAPH_CHECK(igraph_vector_update(res,b));
  }

  if (!cs_ltsolve(L->cs, VECTOR(*res))) {
    IGRAPH_ERROR("Cannot perform lower triangular solve", IGRAPH_FAILURE);
  }
  
  return 0;
}

/**
 * \function igraph_sparsemat_usolve
 * Solve an upper-triangular linear system
 * 
 * Solves the Ux=b upper triangular system.
 * \param U The input matrix, in column-compressed format.
 * \param b The right hand side of the linear system.
 * \param res An initialized vector, the result is stored here.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_usolve(const igraph_sparsemat_t *U,
			    const igraph_vector_t *b,
			    igraph_vector_t *res) {

  if (U->cs->m != U->cs->n) {
    IGRAPH_ERROR("Cannot perform upper triangular solve", IGRAPH_NONSQUARE);
  }

  if (res != b) {
    IGRAPH_CHECK(igraph_vector_update(res, b));
  }

  if (! cs_usolve(U->cs, VECTOR(*res))) {
    IGRAPH_ERROR("Cannot perform upper triangular solve", IGRAPH_FAILURE);
  }
  
  return 0;
}

/**
 * \function igraph_sparsemat_utsolve
 * Solve a lower-triangular linear system
 * 
 * This is the same as \ref igraph_sparsemat_usolve(), but U'x=b is
 * solved, where the apostrophe denotes the transpose.
 * \param U The input matrix, in column-compressed format.
 * \param b The right hand side of the linear system.
 * \param res An initialized vector, the result is stored here.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_utsolve(const igraph_sparsemat_t *U,
			     const igraph_vector_t *b,
			     igraph_vector_t *res) {
  
  if (U->cs->m != U->cs->n) {
    IGRAPH_ERROR("Cannot perform transposed upper triangular solve",
		 IGRAPH_NONSQUARE);
  }

  if (res != b) { 
    IGRAPH_CHECK(igraph_vector_update(res,b));
  }

  if (!cs_utsolve(U->cs, VECTOR(*res))) {
    IGRAPH_ERROR("Cannot perform transposed upper triangular solve", 
		 IGRAPH_FAILURE);
  }
  
  return 0;
}

/** 
 * \function igraph_sparsemat_cholsol
 * Solve a symmetric linear system via Cholesky decomposition
 * 
 * Solve Ax=b, where A is a symmetric positive definite matrix.
 * \param A The input matrix, in column-compressed format.
 * \param v The right hand side.
 * \param res An initialized vector, the result is stored here.
 * \param order An integer giving the ordering method to use for the
 *    factorization. Zero is the natural ordering; if it is one, then
 *    the fill-reducing minimum-degree ordering of A+A' is used.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_cholsol(const igraph_sparsemat_t *A,
			     const igraph_vector_t *b,
			     igraph_vector_t *res, 
			     int order) {
  
  if (A->cs->m != A->cs->n) {
    IGRAPH_ERROR("Cannot perform sparse symmetric solve",
		 IGRAPH_NONSQUARE);
  }

  if (res != b) { 
    IGRAPH_CHECK(igraph_vector_update(res,b));
  }

  if (! cs_cholsol(order, A->cs, VECTOR(*res))) {
    IGRAPH_ERROR("Cannot perform sparse symmetric solve", IGRAPH_FAILURE);
  }
  
  return 0;
}

/**
 * \function igraph_sparsemat_lusol
 * Solve a linear system via LU decomposition
 * 
 * Solve Ax=b, via LU factorization of A.
 * \param A The input matrix, in column-compressed format.
 * \param b The right hand side of the equation.
 * \param res An initialized vector, the result is stored here.
 * \param order The ordering method to use, zero means the natural
 *    ordering, one means the fill-reducing minimum-degree ordering of
 *    A+A', two means the ordering of A'*A, after removing the dense
 *    rows from A. Three means the ordering of A'*A.
 * \param tol Real number, the tolerance limit to use for the numeric
 *    LU factorization.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_lusol(const igraph_sparsemat_t *A,
			   const igraph_vector_t *b,
			   igraph_vector_t *res,
			   int order,
			   igraph_real_t tol) {
  
  if (A->cs->m != A->cs->n) {
    IGRAPH_ERROR("Cannot perform LU solve",
		 IGRAPH_NONSQUARE);
  }

  if (res != b) { 
    IGRAPH_CHECK(igraph_vector_update(res,b));
  }

  if (! cs_lusol(order, A->cs, VECTOR(*res), tol)) {
    IGRAPH_ERROR("Cannot perform LU solve", IGRAPH_FAILURE);
  }
  
  return 0;
}

int igraph_i_sparsemat_cc(igraph_t *graph, const igraph_sparsemat_t *A,
			  igraph_bool_t directed) {

  igraph_vector_t edges;
  long int no_of_nodes=A->cs->m;
  long int no_of_edges=A->cs->p[A->cs->n];
  int *p=A->cs->p;
  int *i=A->cs->i;
  long int from=0;
  long int to=0;
  long int e=0;
  
  if (no_of_nodes != A->cs->n) {
    IGRAPH_ERROR("Cannot create graph object", IGRAPH_NONSQUARE);
  }
  
  IGRAPH_VECTOR_INIT_FINALLY(&edges, no_of_edges*2);
  
  while (*p < no_of_edges) {
    while (to < *(p+1)) {
      if (directed || from >= *i) {
	VECTOR(edges)[e++] = from;
	VECTOR(edges)[e++] = (*i);
      }
      to++;
      i++;      
    }
    from++;
    p++;
  }
  igraph_vector_resize(&edges, e);

  IGRAPH_CHECK(igraph_create(graph, &edges, (igraph_integer_t) no_of_nodes,
			     directed));
  igraph_vector_destroy(&edges);
  IGRAPH_FINALLY_CLEAN(1);

  return 0;
}

int igraph_i_sparsemat_triplet(igraph_t *graph, const igraph_sparsemat_t *A,
			       igraph_bool_t directed) {

  igraph_vector_t edges;
  long int no_of_nodes=A->cs->m;
  long int no_of_edges=A->cs->nz;
  int *i=A->cs->p;
  int *j=A->cs->i;
  long int e;
  
  if (no_of_nodes != A->cs->n) {
    IGRAPH_ERROR("Cannot create graph object", IGRAPH_NONSQUARE);
  }
  
  IGRAPH_VECTOR_INIT_FINALLY(&edges, no_of_edges*2);
  
  for (e=0; e<2*no_of_edges; i++, j++) {
    if (directed || *i >= *j) {
      VECTOR(edges)[e++] = (*i);
      VECTOR(edges)[e++] = (*j);
    }
  }
  igraph_vector_resize(&edges, e);
  
  IGRAPH_CHECK(igraph_create(graph, &edges, (igraph_integer_t) no_of_nodes,
			     directed));
  igraph_vector_destroy(&edges);
  IGRAPH_FINALLY_CLEAN(1);

  return 0;
}

/**
 * \function igraph_sparsemat
 * Create an igraph graph from a sparse matrix
 * 
 * One edge is created for each non-zero entry in the matrix. If you
 * have a symmetric matrix, and want to create an undirected graph,
 * then delete the entries in the upper diagonal first, or call \ref
 * igraph_simplify() on the result graph to eliminate the multiple
 * edges.
 * \param graph Pointer to an uninitialized igraph_t object, the
 *    graphs is stored here.
 * \param A The input matrix, in triplet or column-compressed format.
 * \param directed Boolean scalar, whether to create a directed
 *    graph.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat(igraph_t *graph, const igraph_sparsemat_t *A,
		     igraph_bool_t directed) {
  
  if (A->cs->nz < 0) {
    return(igraph_i_sparsemat_cc(graph, A, directed));
  } else {
    return(igraph_i_sparsemat_triplet(graph, A, directed));
  }
}

int igraph_i_weighted_sparsemat_cc(const igraph_sparsemat_t *A, 
				   igraph_bool_t directed, const char *attr, 
				   igraph_bool_t loops, 
				   igraph_vector_t *edges,
				   igraph_vector_t *weights) {

  long int no_of_edges=A->cs->p[A->cs->n];
  int *p=A->cs->p;
  int *i=A->cs->i;
  igraph_real_t *x=A->cs->x;
  long int from=0;  
  long int to=0;
  long int e=0, w=0;

  IGRAPH_UNUSED(attr);

  igraph_vector_resize(edges, no_of_edges*2);
  igraph_vector_resize(weights, no_of_edges);

  while (*p < no_of_edges) {
    while (to < *(p+1)) {
      if ( (loops || from != *i) && (directed || from >= *i) && *x != 0) {
	VECTOR(*edges)[e++] = (*i);
	VECTOR(*edges)[e++] = from;
	VECTOR(*weights)[w++] = (*x);
      }
      to++;
      i++;
      x++;
    }
    from++;
    p++;
  }

  igraph_vector_resize(edges, e);
  igraph_vector_resize(weights, w);
  
  return 0;
}

int igraph_i_weighted_sparsemat_triplet(const igraph_sparsemat_t *A, 
					igraph_bool_t directed,
					const char *attr, 
					igraph_bool_t loops, 
					igraph_vector_t *edges,
					igraph_vector_t *weights) {

  IGRAPH_UNUSED(A); IGRAPH_UNUSED(directed); IGRAPH_UNUSED(attr);
  IGRAPH_UNUSED(loops); IGRAPH_UNUSED(edges); IGRAPH_UNUSED(weights);

  /* TODO */
  IGRAPH_ERROR("Triplet matrices are not implemented", 
	       IGRAPH_UNIMPLEMENTED);
  return 0;
}

int igraph_weighted_sparsemat(igraph_t *graph, const igraph_sparsemat_t *A, 
			      igraph_bool_t directed, const char *attr, 
			      igraph_bool_t loops) {

  igraph_vector_t edges, weights;
  int pot_edges= A->cs->nz < 0 ? A->cs->p[A->cs->n] : A->cs->nz;
  const char* default_attr = "weight";
  igraph_vector_ptr_t attr_vec;
  igraph_attribute_record_t attr_rec;
  long int no_of_nodes=A->cs->m;

  if (no_of_nodes != A->cs->n) {
    IGRAPH_ERROR("Cannot create graph object", IGRAPH_NONSQUARE);
  }
  
  IGRAPH_VECTOR_INIT_FINALLY(&edges, pot_edges*2);
  IGRAPH_VECTOR_INIT_FINALLY(&weights, pot_edges);
  IGRAPH_VECTOR_PTR_INIT_FINALLY(&attr_vec, 1);

  if (A->cs->nz < 0) {
    IGRAPH_CHECK(igraph_i_weighted_sparsemat_cc(A, directed, attr, loops, 
					       &edges, &weights));
  } else {
    IGRAPH_CHECK(igraph_i_weighted_sparsemat_triplet(A, directed, attr, 
						     loops, &edges, 
						     &weights));
  }

  /* Prepare attribute record */
  attr_rec.name = attr ? attr : default_attr;
  attr_rec.type = IGRAPH_ATTRIBUTE_NUMERIC;
  attr_rec.value = &weights;
  VECTOR(attr_vec)[0] = &attr_rec;

  /* Create graph */
  IGRAPH_CHECK(igraph_empty(graph, (igraph_integer_t) no_of_nodes, directed));
  IGRAPH_FINALLY(igraph_destroy, graph);
  if (igraph_vector_size(&edges)>0) {
    IGRAPH_CHECK(igraph_add_edges(graph, &edges, &attr_vec));
  }
  IGRAPH_FINALLY_CLEAN(1);

  /* Cleanup */
  igraph_vector_destroy(&edges);
  igraph_vector_destroy(&weights);
  igraph_vector_ptr_destroy(&attr_vec);
  IGRAPH_FINALLY_CLEAN(3);
  
  return 0;
}

/**
 * \function igraph_get_sparsemat
 * Convert an igraph graph to a sparse matrix
 * 
 * If the graph is undirected, then a symmetric matrix is created.
 * \param graph The input graph.
 * \param res Pointer to an uninitialized sparse matrix. The result
 *    will be stored here.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_get_sparsemat(const igraph_t *graph, igraph_sparsemat_t *res) {
  
  long int no_of_nodes=igraph_vcount(graph);
  long int no_of_edges=igraph_ecount(graph);
  igraph_bool_t directed=igraph_is_directed(graph);
  long int nzmax= directed ? no_of_edges : no_of_edges*2;
  long int i;  
  
  IGRAPH_CHECK(igraph_sparsemat_init(res, (igraph_integer_t) no_of_nodes, 
				     (igraph_integer_t) no_of_nodes, 
				     (igraph_integer_t) nzmax));
  
  for (i=0; i<no_of_edges; i++) {
    long int from=IGRAPH_FROM(graph, i);
    long int to=IGRAPH_TO(graph, i);
    IGRAPH_CHECK(igraph_sparsemat_entry(res, (int) from, (int) to, 1.0));
    if (!directed && from != to) {
      IGRAPH_CHECK(igraph_sparsemat_entry(res, (int) to, (int) from, 1.0));
    }
  }    
  
  return 0;
}

#define CHECK(x) if ((x)<0) { IGRAPH_ERROR("Cannot write to file", IGRAPH_EFILE); }

/**
 * \function igraph_sparsemat_print
 * Print a sparse matrix to a file
 * 
 * Only the non-zero entries are printed. This function serves more as
 * a debugging utility, as currently there is no function that could
 * read back the printed matrix from the file.
 * \param A The input matrix, triplet or column-compressed format.
 * \param outstream The stream to print it to.
 * \return Error code.
 * 
 * Time complexity: O(nz) for triplet matrices, O(n+nz) for
 * column-compressed matrices. nz is the number of non-zero elements,
 * n is the number columns in the matrix.
 */

int igraph_sparsemat_print(const igraph_sparsemat_t *A,
			   FILE *outstream) {
  
  if (A->cs->nz < 0) {
    /* CC */
    int j, p;
    for (j=0; j<A->cs->n; j++) {
      CHECK(fprintf(outstream, "col %i: locations %i to %i\n", 
		    j, A->cs->p[j], A->cs->p[j+1]-1));
      for (p=A->cs->p[j]; p < A->cs->p[j+1]; p++) {
	CHECK(fprintf(outstream, "%i : %g\n", A->cs->i[p], A->cs->x[p]));
      }
    }
  } else {
    /* Triplet */
    int p;
    for (p=0; p<A->cs->nz; p++) {
      CHECK(fprintf(outstream, "%i %i : %g\n", 
		    A->cs->i[p], A->cs->p[p], A->cs->x[p]));
    }
  }
  
  return 0;
}

#undef CHECK

int igraph_i_sparsemat_eye_triplet(igraph_sparsemat_t *A, int n, int nzmax, 
				   igraph_real_t value) {
  long int i;

  IGRAPH_CHECK(igraph_sparsemat_init(A, n, n, nzmax));
  
  for (i=0; i<n; i++) {
    igraph_sparsemat_entry(A, (int) i, (int) i, value);
  }
  
  return 0;
}

int igraph_i_sparsemat_eye_cc(igraph_sparsemat_t *A, int n,
			      igraph_real_t value) {
  long int i;

  if (! (A->cs = cs_spalloc(n, n, n, /*values=*/ 1, /*triplet=*/ 0)) ) {
    IGRAPH_ERROR("Cannot create eye sparse matrix", IGRAPH_FAILURE);
  }
  
  for (i=0; i<n; i++) {
    A->cs->p [i] = (int) i;
    A->cs->i [i] = (int) i;
    A->cs->x [i] = value;
  }
  A->cs->p [n] = n;
  
  return 0;
}

/**
 * \function igraph_sparsemat_eye
 * Create a sparse identity matrix
 * 
 * \param A An uninitialized sparse matrix, the result is stored
 *   here.
 * \param n The number of rows and number of columns in the matrix.
 * \param nzmax The maximum number of non-zero elements, this
 *   essentially gives the amount of memory that will be allocated for
 *   matrix elements.
 * \param value The value to store in the diagonal.
 * \param compress Whether to create a column-compressed matrix. If
 *   false, then a triplet matrix is created.
 * \return Error code.
 * 
 * Time complexity: O(n).
 */

int igraph_sparsemat_eye(igraph_sparsemat_t *A, int n, int nzmax,
			 igraph_real_t value,
			 igraph_bool_t compress) {
  if (compress) {
    return(igraph_i_sparsemat_eye_cc(A, n, value));
  } else {
    return(igraph_i_sparsemat_eye_triplet(A, n, nzmax, value));
  }
}

int igraph_i_sparsemat_diag_triplet(igraph_sparsemat_t *A, int nzmax,
				    const igraph_vector_t *values) {

  int i, n=(int) igraph_vector_size(values);

  IGRAPH_CHECK(igraph_sparsemat_init(A, n, n, nzmax));
  
  for (i=0; i<n; i++) {
    igraph_sparsemat_entry(A, i, i, VECTOR(*values)[i]);
  }
  
  return 0;
  
}

int igraph_i_sparsemat_diag_cc(igraph_sparsemat_t *A,
			       const igraph_vector_t *values) {

  int i, n=(int) igraph_vector_size(values);

  if (! (A->cs = cs_spalloc(n, n, n, /*values=*/ 1, /*triplet=*/ 0)) ) {
    IGRAPH_ERROR("Cannot create eye sparse matrix", IGRAPH_FAILURE);
  }
  
  for (i=0; i<n; i++) {
    A->cs->p [i] = i;
    A->cs->i [i] = i;
    A->cs->x [i] = VECTOR(*values)[i];
  }
  A->cs->p [n] = n;
  
  return 0;

}

/**
 * \function igraph_sparsemat_diag
 * Create a sparse diagonal matrix
 * 
 * \param A An uninitialized sparse matrix, the result is stored
 *    here.
 * \param nzmax The maximum number of non-zero elements, this
 *   essentially gives the amount of memory that will be allocated for
 *   matrix elements.
 * \param values The values to store in the diagonal, the size of the
 *    matrix defined by the length of this vector.
 * \param compress Whether to create a column-compressed matrix. If
 *   false, then a triplet matrix is created.
 * \return Error code.
 * 
 * Time complexity: O(n), the length of the diagonal vector.
 */

int igraph_sparsemat_diag(igraph_sparsemat_t *A, int nzmax,
			  const igraph_vector_t *values,
			  igraph_bool_t compress) {

  if (compress) {
    return(igraph_i_sparsemat_diag_cc(A, values));
  } else {
    return(igraph_i_sparsemat_diag_triplet(A, nzmax, values));
  }
}

int igraph_i_sparsemat_arpack_multiply(igraph_real_t *to, 
				       const igraph_real_t *from,
				       int n,
				       void *extra) {
  igraph_sparsemat_t *A=extra;
  igraph_vector_t vto, vfrom;
  igraph_vector_view(&vto, to, n);
  igraph_vector_view(&vfrom, from, n);
  igraph_vector_null(&vto);
  IGRAPH_CHECK(igraph_sparsemat_gaxpy(A, &vfrom, &vto));
  return 0;
}

typedef struct igraph_i_sparsemat_arpack_rssolve_data_t {
  igraph_sparsemat_symbolic_t *dis;
  igraph_sparsemat_numeric_t *din;
  igraph_real_t tol;
  igraph_sparsemat_solve_t method;
} igraph_i_sparsemat_arpack_rssolve_data_t;

int igraph_i_sparsemat_arpack_solve(igraph_real_t *to, 
				    const igraph_real_t *from,
				    int n,
				    void *extra) {

  igraph_i_sparsemat_arpack_rssolve_data_t *data=extra;
  igraph_vector_t vfrom, vto;

  igraph_vector_view(&vfrom, from, n);
  igraph_vector_view(&vto, to, n);

  if (data->method == IGRAPH_SPARSEMAT_SOLVE_LU) {
    IGRAPH_CHECK(igraph_sparsemat_luresol(data->dis, data->din, &vfrom, 
					  &vto));
  } else if (data->method == IGRAPH_SPARSEMAT_SOLVE_QR) {
    IGRAPH_CHECK(igraph_sparsemat_qrresol(data->dis, data->din, &vfrom, 
					  &vto));

  }

  return 0;
}

/**
 * \function igraph_sparsemat_arpack_rssolve
 * Eigenvalues and eigenvectors of a symmetric sparse matrix via ARPACK
 * 
 * \param The input matrix, must be column-compressed.
 * \param options It is passed to \ref igraph_arpack_rssolve(). See
 *    \ref igraph_arpack_options_t for the details. If \c mode is 1,
 *    then ARPACK uses regular mode, if \c mode is 3, then shift and
 *    invert mode is used and the \c sigma structure member defines
 *    the shift.
 * \param storage Storage for ARPACK. See \ref
 *    igraph_arpack_rssolve() and \ref igraph_arpack_storage_t for
 *    details.
 * \param values An initialized vector or a null pointer, the
 *    eigenvalues are stored here.
 * \param vectors An initialised matrix, or a null pointer, the
 *    eigenvectors are stored here, in the columns.
 * \param solvemethod The method to solve the linear system, if \c
 *    mode is 3, i.e. the shift and invert mode is used. 
 *    Possible values:
 *    \clist 
 *      \cli IGRAPH_SPARSEMAT_SOLVE_LU
 *           The linear system is solved using LU decomposition.
 *      \cli IGRAPH_SPARSEMAT_SOLVE_QR
 *           The linear system is solved using QR decomposition.
 *    \endclist
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_arpack_rssolve(const igraph_sparsemat_t *A,
				    igraph_arpack_options_t *options,
				    igraph_arpack_storage_t *storage,
				    igraph_vector_t *values,
				    igraph_matrix_t *vectors,
				    igraph_sparsemat_solve_t solvemethod) {

  int n=(int) igraph_sparsemat_nrow(A);

  if (n != igraph_sparsemat_ncol(A)) {
    IGRAPH_ERROR("Non-square matrix for ARPACK", IGRAPH_NONSQUARE);
  }

  options->n=n;

  if (options->mode==1) {
    IGRAPH_CHECK(igraph_arpack_rssolve(igraph_i_sparsemat_arpack_multiply,
				       (void*) A, options, storage, 
				       values, vectors));
  } else if (options->mode==3) {
    igraph_real_t sigma=options->sigma;
    igraph_sparsemat_t OP, eye;
    igraph_sparsemat_symbolic_t symb;
    igraph_sparsemat_numeric_t num;
    igraph_i_sparsemat_arpack_rssolve_data_t data;
    /*-----------------------------------*/
    /* We need to factor the (A-sigma*I) */
    /*-----------------------------------*/
    
    /* Create (A-sigma*I) */
    IGRAPH_CHECK(igraph_sparsemat_eye(&eye, /*n=*/ n, /*nzmax=*/ n, 
				      /*value=*/ -sigma, /*compress=*/ 1));
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &eye);
    IGRAPH_CHECK(igraph_sparsemat_add(/*A=*/ A, /*B=*/ &eye, /*alpha=*/ 1.0,
				      /*beta=*/ 1.0, /*res=*/ &OP));
    igraph_sparsemat_destroy(&eye);
    IGRAPH_FINALLY_CLEAN(1);
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &OP);
        
    if (solvemethod==IGRAPH_SPARSEMAT_SOLVE_LU) {
      /* Symbolic analysis */
      IGRAPH_CHECK(igraph_sparsemat_symblu(/*order=*/ 0, &OP, &symb));
      IGRAPH_FINALLY(igraph_sparsemat_symbolic_destroy, &symb);
      /* Numeric LU factorization */
      IGRAPH_CHECK(igraph_sparsemat_lu(&OP, &symb, &num, /*tol=*/ 0));
      IGRAPH_FINALLY(igraph_sparsemat_numeric_destroy, &num);
    } else if (solvemethod==IGRAPH_SPARSEMAT_SOLVE_QR) {
      /* Symbolic analysis */
      IGRAPH_CHECK(igraph_sparsemat_symbqr(/*order=*/ 0, &OP, &symb));
      IGRAPH_FINALLY(igraph_sparsemat_symbolic_destroy, &symb);
      /* Numeric QR factorization */
      IGRAPH_CHECK(igraph_sparsemat_qr(&OP, &symb, &num));
      IGRAPH_FINALLY(igraph_sparsemat_numeric_destroy, &num);
    }

    data.dis=&symb;
    data.din=&num;
    data.tol=options->tol;
    data.method=solvemethod;
    IGRAPH_CHECK(igraph_arpack_rssolve(igraph_i_sparsemat_arpack_solve,
				       (void*) &data, options, storage, 
				       values, vectors));
   
    igraph_sparsemat_numeric_destroy(&num);
    igraph_sparsemat_symbolic_destroy(&symb);
    igraph_sparsemat_destroy(&OP);
    IGRAPH_FINALLY_CLEAN(3);
  }

  return 0;
}

/**
 * \function igraph_sparsemat_arpack_rnsolve
 * Eigenvalues and eigenvectors of a nonsymmetric sparse matrix via ARPACK
 * 
 * Eigenvalues and/or eigenvectors of a nonsymmetric sparse matrix.
 * \param A The input matrix, in column-compressed mode.
 * \param options ARPACK options, it is passed to \ref
 *    igraph_arpack_rnsolve(). See also \ref igraph_arpack_options_t
 *    for details.
 * \param storage Storage for ARPACK, this is passed to \ref
 *    igraph_arpack_rnsolve(). See \ref igraph_arpack_storage_t for
 *    details.
 * \param values An initialized matrix, or a null pointer. If not a
 *    null pointer, then the eigenvalues are stored here, the first
 *    column is the real part, the second column is the imaginary
 *    part.
 * \param vectors An initialized matrix, or a null pointer. If not a
 *    null pointer, then the eigenvectors are stored here, please see 
 *    \ref igraph_arpack_rnsolve() for the format.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_arpack_rnsolve(const igraph_sparsemat_t *A,
				    igraph_arpack_options_t *options,
				    igraph_arpack_storage_t *storage,
				    igraph_matrix_t *values, 
				    igraph_matrix_t *vectors) {

  int n=(int) igraph_sparsemat_nrow(A);

  if (n != igraph_sparsemat_ncol(A)) {
    IGRAPH_ERROR("Non-square matrix for ARPACK", IGRAPH_NONSQUARE);
  }

  options->n=n;

  return igraph_arpack_rnsolve(igraph_i_sparsemat_arpack_multiply,
			       (void*) A, options, storage, 
			       values, vectors);
}

/**
 * \function igraph_sparsemat_symbqr
 * Symbolic QR decomposition
 * 
 * QR decomposition of sparse matrices involves two steps, the first
 * is calling this function, and then \ref
 * igraph_sparsemat_qr().
 * \param order The ordering to use: 0 means natural ordering, 1 means
 *   minimum degree ordering of A+A', 2 is minimum degree ordering of
 *   A'A after removing the dense rows from A, and 3 is the minimum
 *   degree ordering of A'A. 
 * \param A The input matrix, in column-compressed format.
 * \param dis The result of the symbolic analysis is stored here. Once
 *    not needed anymore, it must be destroyed by calling \ref
 *    igraph_sparsemat_symbolic_destroy().
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_symbqr(long int order, const igraph_sparsemat_t *A,
			    igraph_sparsemat_symbolic_t *dis) {

  dis->symbolic = cs_sqr((int) order, A->cs, /*qr=*/ 1);
  if (!dis->symbolic) {
    IGRAPH_ERROR("Cannot do symbolic QR decomposition", IGRAPH_FAILURE);
  }
  
  return 0;
}

/**
 * \function igraph_sparsemat_symblu
 * Symbolic LU decomposition
 * 
 * LU decomposition of sparse matrices involves two steps, the first
 * is calling this function, and then \ref igraph_sparsemat_lu().
 * \param order The ordering to use: 0 means natural ordering, 1 means
 *   minimum degree ordering of A+A', 2 is minimum degree ordering of
 *   A'A after removing the dense rows from A, and 3 is the minimum
 *   degree ordering of A'A. 
 * \param A The input matrix, in column-compressed format.
 * \param dis The result of the symbolic analysis is stored here. Once
 *    not needed anymore, it must be destroyed by calling \ref
 *    igraph_sparsemat_symbolic_destroy().
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_symblu(long int order, const igraph_sparsemat_t *A,
			    igraph_sparsemat_symbolic_t *dis) {
  
  dis->symbolic = cs_sqr((int) order, A->cs, /*qr=*/ 0);
  if (!dis->symbolic) {
    IGRAPH_ERROR("Cannot do symbolic LU decomposition", IGRAPH_FAILURE);
  }
  
  return 0;
}

/**
 * \function igraph_sparsemat_lu
 * LU decomposition of a sparse matrix
 * 
 * Performs numeric sparse LU decomposition of a matrix.
 * \param A The input matrix, in column-compressed format.
 * \param dis The symbolic analysis for LU decomposition, coming from
 *    a call to the \ref igraph_sparsemat_symblu() function.
 * \param din The numeric decomposition, the result is stored here. It
 *    can be used to solve linear systems with changing right hand
 *    side vectors, by calling \ref igraph_sparsemat_luresol(). Once
 *    not needed any more, it must be destroyed by calling \ref
 *    igraph_sparsemat_symbolic_destroy() on it.
 * \param tol The tolerance for the numeric LU decomposition.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_lu(const igraph_sparsemat_t *A, 
			const igraph_sparsemat_symbolic_t *dis, 
			igraph_sparsemat_numeric_t *din, double tol) {
  din->numeric=cs_lu(A->cs, dis->symbolic, tol);
  if (!din->numeric) {
    IGRAPH_ERROR("Cannot do LU decomposition", IGRAPH_FAILURE);
  }
  return 0;
}

/**
 * \function igraph_sparsemat_qr
 * QR decomposition of a sparse matrix
 * 
 * Numeric QR decomposition of a sparse matrix.
 * \param A The input matrix, in column-compressed format.
 * \param dis The result of the symbolic QR analysis, from the
 *    function \ref igraph_sparsemat_symbqr().
 * \param din The result of the decomposition is stored here, it can
 *    be used to solve many linear systems with the same coefficient
 *    matrix and changing right hand sides, using the \ref
 *    igraph_sparsemat_qrresol() function. Once not needed any more,
 *    one should call \ref igraph_sparsemat_numeric_destroy() on it to
 *    free the allocated memory.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_qr(const igraph_sparsemat_t *A,
			const igraph_sparsemat_symbolic_t *dis,
			igraph_sparsemat_numeric_t *din) {
  din->numeric=cs_qr(A->cs, dis->symbolic);
  if (!din->numeric) {
    IGRAPH_ERROR("Cannot do QR decomposition", IGRAPH_FAILURE);
  }
  return 0;
}

/**
 * \function igraph_sparsemat_luresol
 * Solve linear system using a precomputed LU decomposition
 * 
 * Uses the LU decomposition of a matrix to solve linear systems.
 * \param dis The symbolic analysis of the coefficient matrix, the
 *    result of \ref igraph_sparsemat_symblu().
 * \param din The LU decomposition, the result of a call to \ref
 *    igraph_sparsemat_lu(). 
 * \param b A vector that defines the right hand side of the linear
 *    equation system.
 * \param res An initialized vector, the solution of the linear system
 *    is stored here.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_luresol(const igraph_sparsemat_symbolic_t *dis,
			     const igraph_sparsemat_numeric_t *din, 
			     const igraph_vector_t *b,
			     igraph_vector_t *res) {
  int n=din->numeric->L->n;
  igraph_real_t *workspace;

  if (res != b) {
    IGRAPH_CHECK(igraph_vector_update(res, b));
  }
  
  workspace=igraph_Calloc(n, igraph_real_t);
  if (!workspace) { 
    IGRAPH_ERROR("Cannot LU (re)solve sparse matrix", IGRAPH_ENOMEM);
  }
  IGRAPH_FINALLY(igraph_free, workspace);

  if (!cs_ipvec(din->numeric->pinv, VECTOR(*res), workspace, n)) {
    IGRAPH_ERROR("Cannot LU (re)solve sparse matrix", IGRAPH_FAILURE);
  }
  if (!cs_lsolve(din->numeric->L, workspace)) {
    IGRAPH_ERROR("Cannot LU (re)solve sparse matrix", IGRAPH_FAILURE);
  }
  if (!cs_usolve(din->numeric->U, workspace)) {
    IGRAPH_ERROR("Cannot LU (re)solve sparse matrix", IGRAPH_FAILURE);
  }
  if (!cs_ipvec(dis->symbolic->q, workspace, VECTOR(*res), n)) {
    IGRAPH_ERROR("Cannot LU (re)solve sparse matrix", IGRAPH_FAILURE);
  }
    
  igraph_Free(workspace);
  IGRAPH_FINALLY_CLEAN(1);

  return 0;
}

/**
 * \function igraph_sparsemat_qrresol
 * Solve a linear system using a precomputed QR decomposition
 * 
 * Solves a linear system using a QR decomposition of its coefficient
 * matrix.
 * \param dis Symbolic analysis of the coefficient matrix, the result
 *    of \ref igraph_sparsemat_symbqr().
 * \param din The QR decomposition of the coefficient matrix, the
 *    result of \ref igraph_sparsemat_qr().
 * \param b Vector, giving the right hand side of the linear equation
 *    system. 
 * \param res An initialized vector, the solution is stored here. It
 *    is resized as needed.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_qrresol(const igraph_sparsemat_symbolic_t *dis,
			     const igraph_sparsemat_numeric_t *din, 
			     const igraph_vector_t *b,
			     igraph_vector_t *res) {
  int n=din->numeric->L->n;
  igraph_real_t *workspace;
  int k;

  if (res != b) {
    IGRAPH_CHECK(igraph_vector_update(res, b));
  }

  workspace=igraph_Calloc(dis->symbolic ? dis->symbolic->m2 : 1, 
			  igraph_real_t);
  if (!workspace) { 
    IGRAPH_ERROR("Cannot QR (re)solve sparse matrix", IGRAPH_FAILURE);
  }
  IGRAPH_FINALLY(igraph_free, workspace);

  if (!cs_ipvec(dis->symbolic->pinv, VECTOR(*res), workspace, n)) {
    IGRAPH_ERROR("Cannot QR (re)solve sparse matrix", IGRAPH_FAILURE);
  }
  for (k=0; k<n; k++) {
    if (!cs_happly(din->numeric->L, k, din->numeric->B[k], workspace)) {
      IGRAPH_ERROR("Cannot QR (re)solve sparse matrix", IGRAPH_FAILURE);
    }
  }
  if (!cs_usolve(din->numeric->U, workspace)) {
    IGRAPH_ERROR("Cannot QR (re)solve sparse matrix", IGRAPH_FAILURE);
  }
  if (!cs_ipvec(dis->symbolic->q, workspace, VECTOR(*res), n)) {
    IGRAPH_ERROR("Cannot QR (re)solve sparse matrix", IGRAPH_FAILURE);
  }

  igraph_Free(workspace);
  IGRAPH_FINALLY_CLEAN(1);
  
  return 0;
}

/**
 * \function igraph_sparsemat_symbolic_destroy
 * Deallocate memory for a symbolic decomposition
 * 
 * Frees the memory allocated by \ref igraph_sparsemat_symbqr() or
 * \ref igraph_sparsemat_symblu().
 * \param dis The symbolic analysis.
 * 
 * Time complexity: O(1).
 */

void igraph_sparsemat_symbolic_destroy(igraph_sparsemat_symbolic_t *dis) {
  cs_sfree(dis->symbolic);
  dis->symbolic=0;
}

/**
 * \function igraph_sparsemat_numeric_destroy
 * Deallocate memory for a numeric decomposition
 * 
 * Frees the memoty allocated by \ref igraph_sparsemat_qr() or \ref
 * igraph_sparsemat_lu().
 * \param din The LU or QR decomposition.
 * 
 * Time complexity: O(1).
 */

void igraph_sparsemat_numeric_destroy(igraph_sparsemat_numeric_t *din) {
  cs_nfree(din->numeric);
  din->numeric=0;
}

/**
 * \function igraph_matrix_as_sparsemat
 * Convert a dense matrix to a sparse matrix
 * 
 * \param res An uninitialized sparse matrix, the result is stored
 *    here.
 * \param mat The dense input matrix.
 * \param tol Real scalar, the tolerance. Values closer than \p tol to
 *    zero are considered as zero, and will not be included in the
 *    sparse matrix.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements in the dense
 * matrix.
 */

int igraph_matrix_as_sparsemat(igraph_sparsemat_t *res,
			       const igraph_matrix_t *mat,
			       igraph_real_t tol) {
  int nrow=(int) igraph_matrix_nrow(mat);
  int ncol=(int) igraph_matrix_ncol(mat);
  int i, j, nzmax=0;

  for (i=0; i<nrow; i++) {
    for (j=0; j<ncol; j++) {
      if (fabs(MATRIX(*mat, i, j)) > tol) { nzmax++; }
    }
  }
  
  IGRAPH_CHECK(igraph_sparsemat_init(res, nrow, ncol, nzmax));
  
  for (i=0; i<nrow; i++) {
    for (j=0; j<ncol; j++) {
      if (fabs(MATRIX(*mat, i, j)) > tol) {
	IGRAPH_CHECK(igraph_sparsemat_entry(res, i, j, MATRIX(*mat, i, j)));
      }
    }
  }

  return 0;
}

int igraph_i_sparsemat_as_matrix_cc(igraph_matrix_t *res,
				    const igraph_sparsemat_t *spmat) {

  int nrow=(int) igraph_sparsemat_nrow(spmat);
  int ncol=(int) igraph_sparsemat_ncol(spmat);
  int *p=spmat->cs->p;
  int *i=spmat->cs->i;
  igraph_real_t *x=spmat->cs->x;
  int nzmax=spmat->cs->nzmax;
  int from=0, to=0;

  IGRAPH_CHECK(igraph_matrix_resize(res, nrow, ncol));
  igraph_matrix_null(res);
  
  while (*p < nzmax) {
    while (to < *(p+1)) {
      MATRIX(*res, *i, from) += *x;
      to++;
      i++;
      x++;
    }
    from++;
    p++;
  }

  return 0;
}

int igraph_i_sparsemat_as_matrix_triplet(igraph_matrix_t *res,
					 const igraph_sparsemat_t *spmat) {
  int nrow=(int) igraph_sparsemat_nrow(spmat);
  int ncol=(int) igraph_sparsemat_ncol(spmat);
  int *i=spmat->cs->p;
  int *j=spmat->cs->i;
  igraph_real_t *x=spmat->cs->x;
  int nz=spmat->cs->nz;
  int e;
  
  IGRAPH_CHECK(igraph_matrix_resize(res, nrow, ncol));
  igraph_matrix_null(res);

  for (e=0; e<nz; e++, i++, j++, x++) {
    MATRIX(*res, *j, *i) += *x;
  }

  return 0;
}

/**
 * \function igraph_sparsemat_as_matrix
 * Convert a sparse matrix to a dense matrix
 * 
 * \param res Pointer to an initialized matrix, the result is stored
 *    here. It will be resized to the required size.
 * \param spmat The input sparse matrix, in triplet or
 *    column-compressed format.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements in the dense
 * matrix.
 */

int igraph_sparsemat_as_matrix(igraph_matrix_t *res,
			       const igraph_sparsemat_t *spmat) {
  if (spmat->cs->nz < 0) {
    return(igraph_i_sparsemat_as_matrix_cc(res, spmat));
  } else {
    return(igraph_i_sparsemat_as_matrix_triplet(res, spmat));
  }
}

/**
 * \function igraph_sparsemat_max
 * Maximum of a sparse matrix
 * 
 * \param A The input matrix, column-compressed.
 * \return The maximum in the input matrix, or \c IGRAPH_NEGINFINITY
 *    if the matrix has zero elements.
 * 
 * Time complexity: TODO.
 */

igraph_real_t igraph_sparsemat_max(igraph_sparsemat_t *A) {
  int i, n;
  igraph_real_t *ptr;
  igraph_real_t res;

  IGRAPH_CHECK(igraph_sparsemat_dupl(A));

  ptr=A->cs->x;
  n = A->cs->nz==-1 ? A->cs->p[A->cs->n] : A->cs->nz;
  if (n==0) { return IGRAPH_NEGINFINITY;  } 
  res = *ptr;
  for (i=1; i<n; i++, ptr++) {
    if (*ptr > res) { res=*ptr; }
  }
  return res;
}

/* TODO: CC matrix don't actually need _dupl, 
   because the elements are right beside each other. 
   Same for max and minmax. */

/**
 * \function igraph_sparsemat_min
 * Minimum of a sparse matrix
 *
 * \param A The input matrix, column-compressed.
 * \return The minimum in the input matrix, or \c IGRAPH_POSINFINITY
 *    if the matrix has zero elements.
 * 
 * Time complexity: TODO.
 */

igraph_real_t igraph_sparsemat_min(igraph_sparsemat_t *A) {
  int i, n;
  igraph_real_t *ptr;
  igraph_real_t res;

  IGRAPH_CHECK(igraph_sparsemat_dupl(A));

  ptr=A->cs->x;
  n = A->cs->nz==-1 ? A->cs->p[A->cs->n] : A->cs->nz;
  if (n==0) { return IGRAPH_POSINFINITY;  } 
  res = *ptr;
  for (i=1; i<n; i++, ptr++) {
    if (*ptr < res) { res=*ptr; }
  }
  return res;
}

/**
 * \function igraph_sparsemat_minmax
 * Minimum and maximum of a sparse matrix
 *
 * \param A The input matrix, column-compressed.
 * \param min The minimum in the input matrix is stored here, or \c
 *    IGRAPH_POSINFINITY if the matrix has zero elements.
 * \param max The maximum in the input matrix is stored here, or \c
 *    IGRAPH_NEGINFINITY if the matrix has zero elements.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */


int igraph_sparsemat_minmax(igraph_sparsemat_t *A, 
			    igraph_real_t *min, igraph_real_t *max) {
  int i, n;
  igraph_real_t *ptr;

  IGRAPH_CHECK(igraph_sparsemat_dupl(A));

  ptr=A->cs->x;
  n = A->cs->nz==-1 ? A->cs->p[A->cs->n] : A->cs->nz;
  if (n==0) { 
    *min=IGRAPH_POSINFINITY;
    *max=IGRAPH_NEGINFINITY;
    return 0;
  }
  *min = *max = *ptr;
  for (i=1; i<n; i++, ptr++) {
    if (*ptr > *max) { 
      *max=*ptr; 
    } else if (*ptr < *min) {
      *min=*ptr;
    }
  }
  return 0;
}

/**
 * \function igraph_sparsemat_count_nonzero
 * Count nonzero elements of a sparse matrix
 * 
 * \param A The input matrix, column-compressed.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

long int igraph_sparsemat_count_nonzero(igraph_sparsemat_t *A) {
  int i, n;
  int res=0;
  igraph_real_t *ptr;

  IGRAPH_CHECK(igraph_sparsemat_dupl(A));

  ptr=A->cs->x;
  n = A->cs->nz==-1 ? A->cs->p[A->cs->n] : A->cs->nz;
  if (n==0) { return 0; }
  for (i=0; i<n; i++, ptr++) {
    if (*ptr) { res++; }
  }
  return res;
}

/**
 * \function igraph_sparsemat_count_nonzerotol
 * Count nonzero elements of a sparse matrix, ignoring elements close to zero
 * 
 * Count the number of matrix entries that are closer to zero than \p
 * tol.
 * \param The input matrix, column-compressed.
 * \param Real scalar, the tolerance. 
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

long int igraph_sparsemat_count_nonzerotol(igraph_sparsemat_t *A, 
					   igraph_real_t tol) {
  int i, n;
  int res=0;
  igraph_real_t *ptr;

  IGRAPH_CHECK(igraph_sparsemat_dupl(A));

  ptr=A->cs->x;
  n = A->cs->nz==-1 ? A->cs->p[A->cs->n] : A->cs->nz;
  if (n==0) { return 0; }
  for (i=0; i<n; i++, ptr++) {
    if (*ptr < - tol || *ptr > tol) { res++; }
  }
  return res;
}

int igraph_i_sparsemat_rowsums_triplet(const igraph_sparsemat_t *A,
				       igraph_vector_t *res) {
  int i;
  int *pi=A->cs->i;
  double *px=A->cs->x;
  
  IGRAPH_CHECK(igraph_vector_resize(res, A->cs->m));
  igraph_vector_null(res);
  
  for (i=0; i<A->cs->nz; i++, pi++, px++) {
    VECTOR(*res)[ *pi ] += *px;
  }
  
  return 0;
}

int igraph_i_sparsemat_rowsums_cc(const igraph_sparsemat_t *A,
				  igraph_vector_t *res) {
  int ne=A->cs->p[A->cs->n];
  double *px=A->cs->x;
  int *pi=A->cs->i;

  IGRAPH_CHECK(igraph_vector_resize(res, A->cs->m));
  igraph_vector_null(res);
  
  for (; pi < A->cs->i+ne; pi++, px++) {
    VECTOR(*res)[ *pi ] += *px;
  }

  return 0;
}

/**
 * \function igraph_sparsemat_rowsums
 * Row-wise sums.
 * 
 * \param A The input matrix, in triplet or column-compressed format.
 * \param res An initialized vector, the result is stored here. It
 *    will be resized as needed.
 * \return Error code.
 * 
 * Time complexity: O(nz), the number of non-zero elements.
 */

int igraph_sparsemat_rowsums(const igraph_sparsemat_t *A, 
			     igraph_vector_t *res) {
  if (igraph_sparsemat_is_triplet(A)) {
    return igraph_i_sparsemat_rowsums_triplet(A, res);
  } else {
    return igraph_i_sparsemat_rowsums_cc(A, res);
  }
}

int igraph_i_sparsemat_colsums_triplet(const igraph_sparsemat_t *A,
				       igraph_vector_t *res) {
  int i;
  int *pp=A->cs->p;
  double *px=A->cs->x;
  
  IGRAPH_CHECK(igraph_vector_resize(res, A->cs->n));
  igraph_vector_null(res);
  
  for (i=0; i<A->cs->nz; i++, pp++, px++) {
    VECTOR(*res)[ *pp ] += *px;
  }
  
  return 0;
}

int igraph_i_sparsemat_colsums_cc(const igraph_sparsemat_t *A,
				  igraph_vector_t *res) {
  int n=A->cs->n;
  double *px=A->cs->x;
  int *pp=A->cs->p;
  int *pi=A->cs->i;
  double *pr;
  
  IGRAPH_CHECK(igraph_vector_resize(res, n));
  igraph_vector_null(res);
  pr=VECTOR(*res);
  
  for (; pp < A->cs->p + n; pp++, pr++) {
    for (; pi < A->cs->i + *(pp+1); pi++, px++) {
      *pr += *px;
    }
  }
  return 0;
}

/**
 * \function igraph_sparsemat_colsums
 * Column-wise sums
 * 
 * \param A The input matrix, in triplet or column-compressed format.
 * \param res An initialized vector, the result is stored here. It
 *    will be resized as needed.
 * \return Error code.
 * 
 * Time complexity: O(nz) for triplet matrices, O(nz+n) for
 * column-compressed ones, nz is the number of non-zero elements, n is
 * the number of columns.
 */

int igraph_sparsemat_colsums(const igraph_sparsemat_t *A, 
			     igraph_vector_t *res) {
  if (igraph_sparsemat_is_triplet(A)) {
    return igraph_i_sparsemat_colsums_triplet(A, res);
  } else {
    return igraph_i_sparsemat_colsums_cc(A, res);
  }
}

/**
 * \function igraph_sparsemat_scale
 * Scale a sparse matrix
 * 
 * Multiplies all elements of a sparse matrix, by the given scalar.
 * \param A The input matrix.
 * \param by The scaling factor.
 * \return Error code.
 * 
 * Time complexity: O(nz), the number of non-zero elements in the
 * matrix.
 */

int igraph_sparsemat_scale(igraph_sparsemat_t *A, igraph_real_t by) {
			   
  double *px = A->cs->x;
  int n = A->cs->nz == -1 ? A->cs->p[A->cs->n] : A->cs->nz;
  double *stop=px+n;
  
  for (; px < stop; px++) {
    *px *= by;
  }
  
  return 0;
}

/** 
 * \function igraph_sparsemat_add_rows
 * Add rows to a sparse matrix
 * 
 * The current matrix elements are retained and all elements in the
 * new rows are zero.
 * \param A The input matrix, in triplet or column-compressed format.
 * \param n The number of rows to add.
 * \return Error code.
 * 
 * Time complexity: O(1).
 */

int igraph_sparsemat_add_rows(igraph_sparsemat_t *A, long int n) {
  A->cs->m += n;
  return 0;
}

/**
 * \function igraph_sparsemat_add_cols
 * Add columns to a sparse matrix
 * 
 * The current matrix elements are retained, and all elements in the
 * new columns are zero.
 * \param A The input matrix, in triplet or column-compressed format.
 * \param n The number of columns to add.
 * \return Error code.
 * 
 * Time complexity: TODO.
 */

int igraph_sparsemat_add_cols(igraph_sparsemat_t *A, long int n) {
  if (igraph_sparsemat_is_triplet(A)) {
    A->cs->n += n;
  } else {
    int *newp=realloc(A->cs->p, sizeof(int) * (size_t) (A->cs->n + n + 1));
    int i;
    if (!newp) { 
      IGRAPH_ERROR("Cannot add columns to sparse matrix", IGRAPH_ENOMEM);
    }
    if (newp != A->cs->p) { 
      A->cs->p=newp;
    }
    for (i=A->cs->n+1; i<A->cs->n + n + 1; i++) {
      A->cs->p[i]=A->cs->p[i-1];
    }
    A->cs->n += n;
  }
  return 0;
}

/**
 * \function igraph_sparsemat_resize
 * Resize a sparse matrix
 * 
 * This function resizes a sparse matrix. The resized sparse matrix 
 * will be empty.
 * 
 * \param A The initialized sparse matrix to resize.
 * \param nrow The new number of rows.
 * \param ncol The new number of columns.
 * \param nzmax The new maximum number of elements.
 * \return Error code.
 * 
 * Time complexity: O(nzmax), the maximum number of non-zero elements.
 */

int igraph_sparsemat_resize(igraph_sparsemat_t *A, long int nrow, 
			    long int ncol, int nzmax) {

  if (A->cs->nz < 0) {
    igraph_sparsemat_t tmp;
    IGRAPH_CHECK(igraph_sparsemat_init(&tmp, (int) nrow, (int) ncol, nzmax));
    igraph_sparsemat_destroy(A);
    *A = tmp;
  } else {
    IGRAPH_CHECK(igraph_sparsemat_realloc(A, nzmax));
    A->cs->m = (int) nrow;
    A->cs->n = (int) ncol;
    A->cs->nz = 0;
  }
  return 0;
}

int igraph_sparsemat_nonzero_storage(const igraph_sparsemat_t *A) {
  if (A->cs->nz < 0) {
    return A->cs->p[A->cs->n];
  } else {
    return A->cs->nz;
  }
}

int igraph_sparsemat_getelements(const igraph_sparsemat_t *A, 
				igraph_vector_int_t *i,
				igraph_vector_int_t *j, 
				igraph_vector_t *x) {
  int nz=A->cs->nz;
  if (nz < 0) {
    nz=A->cs->p[A->cs->n];
    IGRAPH_CHECK(igraph_vector_int_resize(i, nz));
    IGRAPH_CHECK(igraph_vector_int_resize(j, A->cs->n+1));
    IGRAPH_CHECK(igraph_vector_resize(x, nz));
    memcpy(VECTOR(*i), A->cs->i, (size_t) nz * sizeof(int));
    memcpy(VECTOR(*j), A->cs->p, (size_t) (A->cs->n+1) * sizeof(int));
    memcpy(VECTOR(*x), A->cs->x, (size_t) nz * sizeof(igraph_real_t));
  } else {
    IGRAPH_CHECK(igraph_vector_int_resize(i, nz));
    IGRAPH_CHECK(igraph_vector_int_resize(j, nz));
    IGRAPH_CHECK(igraph_vector_resize(x, nz));
    memcpy(VECTOR(*i), A->cs->i, (size_t) nz * sizeof(int));
    memcpy(VECTOR(*j), A->cs->p, (size_t) nz * sizeof(int));
    memcpy(VECTOR(*x), A->cs->x, (size_t) nz * sizeof(igraph_real_t));
  }
  return 0;
}

int igraph_sparsemat_scale_rows(igraph_sparsemat_t *A,
				const igraph_vector_t *fact) {
  int *i=A->cs->i;
  igraph_real_t *x=A->cs->x;
  int no_of_edges=A->cs->nz < 0 ? A->cs->p[A->cs->n] : A->cs->nz;
  int e;

  for (e=0; e<no_of_edges; e++, x++, i++) {
    igraph_real_t f=VECTOR(*fact)[*i];
    (*x) *= f;
  }

  return 0;
}

int igraph_i_sparsemat_scale_cols_cc(igraph_sparsemat_t *A,
				     const igraph_vector_t *fact) {
  int *i=A->cs->i;
  igraph_real_t *x=A->cs->x;
  int no_of_edges=A->cs->p[A->cs->n];
  int e;
  int c=0;			/* actual column */
  
  for (e=0; e<no_of_edges; e++, x++, i++) {
    igraph_real_t f;
    while (c < A->cs->n && A->cs->p[c+1] == e) { c++; }
    f=VECTOR(*fact)[c];
    (*x) *= f;
  }

  return 0;
}

int igraph_i_sparsemat_scale_cols_triplet(igraph_sparsemat_t *A,
					  const igraph_vector_t *fact) {
  int *j=A->cs->p;
  igraph_real_t *x=A->cs->x;
  int no_of_edges=A->cs->nz;
  int e;

  for (e=0; e<no_of_edges; e++, x++, j++) {
    igraph_real_t f=VECTOR(*fact)[*j];
    (*x) *= f;
  }

  return 0;
}

int igraph_sparsemat_scale_cols(igraph_sparsemat_t *A,
				const igraph_vector_t *fact) {
  if (A->cs->nz < 0) {
    return igraph_i_sparsemat_scale_cols_cc(A, fact);
  } else {
    return igraph_i_sparsemat_scale_cols_triplet(A, fact);
  }
}

int igraph_sparsemat_multiply_by_dense(const igraph_sparsemat_t *A,
				       const igraph_matrix_t *B,
				       igraph_matrix_t *res) {

  int m=(int) igraph_sparsemat_nrow(A);
  int n=(int) igraph_sparsemat_ncol(A);
  int p=(int) igraph_matrix_ncol(B);
  int i;
  
  if (igraph_matrix_nrow(B) != n) {
    IGRAPH_ERROR("Invalid dimensions in sparse-dense matrix product",
		 IGRAPH_EINVAL);
  }
  
  IGRAPH_CHECK(igraph_matrix_resize(res, m, p));
  igraph_matrix_null(res);

  for (i=0; i<p; i++) {
    if (!(cs_gaxpy(A->cs, &MATRIX(*B, 0, i), &MATRIX(*res, 0, i)))) {
      IGRAPH_ERROR("Cannot perform sparse-dense matrix multiplication",
		   IGRAPH_FAILURE);
    }
  }
  
  return 0;
}

int igraph_sparsemat_dense_multiply(const igraph_matrix_t *A,
				    const igraph_sparsemat_t *B,
				    igraph_matrix_t *res) {
  int m=(int) igraph_matrix_nrow(A);
  int n=(int) igraph_matrix_ncol(A);
  int p=(int) igraph_sparsemat_ncol(B);
  int r, c;
  int *Bp=B->cs->p;
  
  if (igraph_sparsemat_nrow(B) != n) {
    IGRAPH_ERROR("Invalid dimensions in dense-sparse matrix product",
		 IGRAPH_EINVAL);
  }

  if (!igraph_sparsemat_is_cc(B)) {
    IGRAPH_ERROR("Dense-sparse product is only implemented for "
		 "column-compressed sparse matrices", IGRAPH_EINVAL);
  }
  
  IGRAPH_CHECK(igraph_matrix_resize(res, m, p));
  igraph_matrix_null(res);
  
  for (c=0; c<p; c++) {
    for (r=0; r<m; r++) {
      int idx=*Bp;
      while (idx < *(Bp+1)) {
	MATRIX(*res, r, c) += MATRIX(*A, r, B->cs->i[idx]) * B->cs->x[idx];
	idx++;
      }
    }
    Bp++;
  }
  
  return 0;
}

int igraph_i_sparsemat_view(igraph_sparsemat_t *A, int nzmax, int m, int n, 
			    int *p, int *i, double *x, int nz) {

  A->cs = cs_calloc(1, sizeof(cs_di));
  A->cs->nzmax = nzmax;
  A->cs->m = m;
  A->cs->n = n;
  A->cs->p = p;
  A->cs->i = i;
  A->cs->x = x;
  A->cs->nz = nz;
  
  return 0;
}

int igraph_sparsemat_sort(const igraph_sparsemat_t *A, 
			  igraph_sparsemat_t *sorted) {
  
  igraph_sparsemat_t tmp;

  IGRAPH_CHECK(igraph_sparsemat_transpose(A, &tmp, /*values=*/ 1));
  IGRAPH_FINALLY(igraph_sparsemat_destroy, &tmp);
  IGRAPH_CHECK(igraph_sparsemat_transpose(&tmp, sorted, /*values=*/ 1));
  igraph_sparsemat_destroy(&tmp);
  IGRAPH_FINALLY_CLEAN(1);

  return 0;
}

int igraph_sparsemat_getelements_sorted(const igraph_sparsemat_t *A, 
					igraph_vector_int_t *i,
					igraph_vector_int_t *j, 
					igraph_vector_t *x) {
  if (A->cs->nz < 0) {
    igraph_sparsemat_t tmp;
    IGRAPH_CHECK(igraph_sparsemat_sort(A, &tmp));
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &tmp);
    IGRAPH_CHECK(igraph_sparsemat_getelements(&tmp, i, j, x));
    igraph_sparsemat_destroy(&tmp);
    IGRAPH_FINALLY_CLEAN(1);
  } else {
    IGRAPH_CHECK(igraph_sparsemat_getelements(A, i, j, x));
  }

  return 0;
}

int igraph_sparsemat_nzmax(const igraph_sparsemat_t *A) {
  return A->cs->nzmax;
}

int igraph_sparsemat_neg(igraph_sparsemat_t *A) {
  int i, nz=A->cs->nz == -1 ? A->cs->p[A->cs->n] : A->cs->nz;
  igraph_real_t *px=A->cs->x;

  for (i=0; i<nz; i++, px++) { 
    *px = - (*px);
  }
  
  return 0;
}

int igraph_sparsemat_iterator_init(igraph_sparsemat_iterator_t *it,
				   igraph_sparsemat_t *sparsemat) {

  it->mat=sparsemat;
  igraph_sparsemat_iterator_reset(it);
  return 0;
}

int igraph_sparsemat_iterator_reset(igraph_sparsemat_iterator_t *it) {
  it->pos=0;
  if (!igraph_sparsemat_is_triplet(it->mat)) {
    it->col=0;
    while (it->col < it->mat->cs->n && 
	   it->mat->cs->p[it->col+1] == it->pos) { it->col ++; }
  }
  return 0;
}

igraph_bool_t 
igraph_sparsemat_iterator_end(const igraph_sparsemat_iterator_t *it) {
  int nz=it->mat->cs->nz == -1 ? it->mat->cs->p[it->mat->cs->n] : 
    it->mat->cs->nz;
  return it->pos >= nz;
}

int igraph_sparsemat_iterator_row(const igraph_sparsemat_iterator_t *it) {
  return it->mat->cs->i[it->pos];
}

int igraph_sparsemat_iterator_col(const igraph_sparsemat_iterator_t *it) {
  if (igraph_sparsemat_is_triplet(it->mat)) {
    return it->mat->cs->p[it->pos];
  } else {
    return it->col;
  }
}

igraph_real_t 
igraph_sparsemat_iterator_get(const igraph_sparsemat_iterator_t *it) {
  return it->mat->cs->x[it->pos];
}

int igraph_sparsemat_iterator_next(igraph_sparsemat_iterator_t *it) {
  it->pos += 1;
  while (it->col < it->mat->cs->n && 
	 it->mat->cs->p[it->col+1] == it->pos) {
    it->col++;
  }
  return it->pos;
}

int igraph_sparsemat_iterator_idx(const igraph_sparsemat_iterator_t *it) {
  return it->pos;
}
